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ABSTRACT 

A symplectic integrator algorithm suitable for hierarchical triple systems is formulated and 
tested. The positions of the stars are followed in hierarchical Jacobi coordinates, whilst the 
planets are referenced purely to their primary. The algorithm is fast, accurate and easily gen- 
eralised to incorporate collisions. There are five distinct cases - circumtriple orbits, circumbi- 
nary orbits and circumstellar orbits around each of the stars in the hierarchical triple - which 
require a different formulation of the symplectic integration algorithm. 

As an application, a survey of the stability zones for planets in hierarchical triples is 
presented, with the case of a single planet orbiting the inner binary considered in detail. Fits 
to the inner and outer edges of the stability zone are computed. Considering the hierarchical 
triple as two decoupled binary systems, the earlier work of Holman & Wiegert on binaries is 
shown to be applicable to triples, except in the cases of high eccentricities and close or massive 
stars. Application to triple stars with good data in the multiple star catalogue suggests that 
more than 50 % are unable to support circumbinary planets, as the stable zone is non-existent 
or very narrow. 

Key words: celestial mechanics - planetary systems - binaries: general - methods: A^-body 
simulations 



1 INTRODUCTION 

It is well known that many of the stars in the solar neighbourhood exist in multiple systems. As the number of planetary surveys increases, 
planets are regularly being found not only in single star systems, but binaries and triples as well. F or example, recently a hot J upiter has been 



claimed in the triple system HD 188753 (Konacki 2005; but see also Eggenberger et al. 2007), and lDesidera & Barbieril ( 120071) lists 33 binary 
systems and 2 other hierarchical triples known to harbour exoplanets. As the majority of work on planetary dynamics has been for single 
star systems, the dynamics of bodies in these multiple stellar systems is of great interest. At first sight, it might not seem likely to expect 
long term stable planetary systems to exist in binary star systems, let alone triples, but numerical and observational work is starting to show 
otherwise. 

In recent years, much study has been devoted to the stabili ty of planets and planetesimal discs in bi nary star systems. There are several 
investigations of individual systems (e.g. work on 7 Cephei bv D vorak et aljEo03 . Haghighipoui 20061 and lVerrie r & Evans"200d) as wel l 



as substantial amounts of work on more general limits for stability of smaller bodies in these systems. Notablv. Eolman & Wieeert 1 1999h 



approach this problem by using numerical simulation data to empirical fit critical semimajor axes for test particle stability as functions of 
binary mass ratio and eccentricity. These general studies are of great use in the investigation of observed systems and their stability properties, 
giving an effective and fast method of placing limits on stability in the large parameter space created by observational uncertainties. 

The aim of this work is to investigate test particle stability in hierarchical triple star systems, and to see if any similar boundaries can be 
defined for them. To do this, a new method for numerically integrating planets in such systems is presented, following the ideas for binary 
systems presented by Chambers et al. 1 20021) . 



Although there have been empirical studies of the stability of hierarchical star systems themselves, there do not appear to have been 
studies of small b odies in such sy stems. There is a great deal of literature concerning periodic orbits in the general three and four body 



problems (see e.g. I Szebehelvll 1 9671) . but these are almost entirely devoted to considerations of planetary satellites in single star systems, for 
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example satellite transfer in the Sun-Earth-Moon systems. Also, while periodic orbits prove that stable solutions exist in these problems, they 
are of little practical use in determining general stability limits. 

The problem of planetary orbits in triple systems is more complex than for those in binary systems, as there are many different orbital 
configurations possible relative to the three stars. These are considered in Section[2l and a method for classifying them is described in order to 
simplify the discussion in this paper. The derivation of a method to numerically integrate such planets is then given in Section[3] In Section|4] 
is a brief overview of the statistical properties of known triple systems, as a basis for deciding the parameters of the systems used in the 
numerical simulations. The results of numerical investigations into stability properties are then presented in Section[5]for one of the possible 
orbital types. 



2 THE TYPES OF STELLAR AND PLANETARY ORBITS 



The triple systems studied here are all hierarchical in nature. That is, they can all be considered to be a close binary orbited by a distant 
companion - in effect, an inner binary pair and an outer binary pair. Here, the inner binary stars are labelled as A and B, with A being the 
more massive, and the distan t comp anion star as C. There are however three possible ways to pair the three stars into the hierarchy described 
above. Eggleton & Kiseleval i 19951) define the hierarchy by requiring that the instantaneous Keplerian orbits of the inner and outer pair are 



bound, that the outer binary has a longer period than the inner binary and that the ratio of the two periods is the largest of the three possible 
pairings of the stars. This definition is adopted here. 

Although other, non-hierarchical, ty pes of motion are possible for triple stars (see for example Szebehelv 1977 *) and observed (see for 
example the trapezium systems listed in TokovininI 1997 ), they are not considered here. It is an open question whether non-hierarchical 
systems can sustain long-term stable planets. 

Next some system of classifying planetary orbits is needed. The orbital motion of small particles under the influence of two or three 
masses has been studied to some extent through periodic o rbits of the three a nd four body problems, as mentioned in the Introduction. As a 
result of this many families and classification schemes exist. lSzebehelvl j 19671) gives a good review, mentioning for example the (a) to (r) types 
designated for the Copenhagen problem, dependent on the nature and location of the particles motion in both inertial and corotating frames. 
Many other features have been used to desc ribe orbital famili es, for example symmetry (,Brouck&,2004,) or a parameter used to generate 
the orbital family. However, as mentioned bv lSzebehelvl ( Il967h . there is no overall method, and the descriptions would not be applicable to 
non-periodic orbits. 

Orbital types in binary systems have generally been classified into three main groups, those around a single star (circumstellar), those 
around both stars (circumbinary) and those in the middle ground as it were, coorbiting with one of the stars i.e. l ibrating about t he triangular 
Lagrangian points, similar to Trojan asteroids in the Solar System. A convenient labelling of these was given bv lPvorakl Jl984 '), who desig- 
nated planets as either S-orbits (satellite), P-orbits (planet), or L-orbits (librator) for circumstellar, circumbinary and coorbital respectively. 
These ideas can be extended to hierarchical triple systems fairly simply. It is clear that there will be three types of circumstellar motion, one 
for each star, and as such orbits like this can be labelled S(A), S(B) and S(C) for their primary. The circumtriple orbit about the centre of 
mass of all three can be identified as a P orbit and labelled as such in analogy to Dvorak's scheme. Finally, planets which orbit the centre 
of mass of the inner binary are circumbinary but also share characteristics with the satellite orbit, and can be labelled with the combination 
S(AB)-P. These orbital types are shown in Figure [T] along with the binary cases for comparison, and listed in Table[T] A superscript 2 has 
been given to those in binary systems and a 3 for those in triple systems for clarity. 

Although these are the only types of motion studied here, for completeness the coorbital types are also be included. This type of motion 
can occur for both the inner and outer binary, and labels (AB) and (ABC) can be used to designate this. Instead of the single L-orbit type they 
can be broken up into T and H orbits to indicated tadpole (about one of the triangular Lagrangian points) or horseshoe (about both triangular 
points and the L3 collinear point). These are again illustrated in Figure[T]and included in Table[T] 

These orbits are defined to be bound relative to their focus and hierarchical in the same sense as the definition given for the stellar 
system. A particle that orbits outside the extent of the inner binary but has a bound orbit with respect to star A and the binary centre of mass 
would be an S(AB)-P orbit and not an S(A) orbit. 

This method of labelling clearly and concisely designates the exact nature of a planetary orbit, and extends a system already in use for 
binary stars. It can also be easily applied to systems with other levels of hierarchy, for example if star C was replaced with another close 
binary pair C and D additional classes S(D) and S(CD)-P would be possible, and the outer P type could be relabelled P(ABCD). 



3 A SYMPLECTIC INTEGRATOR FOR HIERARCHICAL TRIPLES 

Given the number of simulations that are needed to accurately describe the stability of general triple systems, an accurate and fast method of 
numerically integrating the equations of motion is required. Methods such as Bulirsch-Stoer are reasonably accurate, but very slow. Instead, 
a symplectic method is favoured, as it is fast and shows excellent energy conserv ation. 

A symplectic integration method for planetary systems was introduced by Wisdom & Holmai] 1 199 ih . Here, the Hamiltonian of the 
system is split into a dominant Keplerian part and smaller interaction terms, all analytically integrable on their own. A leapfrog scheme is 
then applied to evolve the system forwards by one timestep. If Hk is the Keplerian Hamiltonian and Hi the smaller interaction terms, then 
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Figure 1. Basic planetary orbit types in multiple star systems as described in Section |2] The top row shows a binary system, in a frame corotating with the 
stars. The bottom row shows a triple system, in a frame corotating with either the inner or outer binary as appropriate. Stellar orbits are shown with dashed 
lines and stars marked with an asteiisk. Pla netary orbits are s hown with a solid line. Note that these plots show real data from numer ical simulations generated 
using a standard Bulirsch-Stoer integrator jPress et alj|2002l) for the tadpole and horseshoe orbits (based on initial conditions from [Murray & Dermotlll200Cll) 
and the integration scheme presented in this paper for the others. 




Figure 2. The coordinate system used for the symplectic integrator. Planets are all taken as relative to their primary, whilst each planetary-stellar subsystem 
barycentre is referred to the centre of mass of the preceding objects. 



to move forward by one step of dt the system is evolved by Hj for Hk for dt and finally by Hj for again. This is a 2nd order 

method. 

There have been two symplectic integrators derived for planets in multiple star systems, by Chambers et al. ( 2002 ) and Beust 1 200 3h 
respectively. IBcusI i 200 3 ^) uses modified Jacobi coordinates for hierarchical systems of any multiplicity. These 'hierarchical Jacobi coordi- 
nates' account for the separate substructures in the hierarchy. Following the example given by Beust ( 2003.), if a system consists of four stars 
in two separate binaries orbiting the systems ba rycentre, then each pai r is assigned Jacobi coordinates within their subsystems, and then Ja- 
cobi coordinates are applied to the two binaries. [Chambers et d] i2002h uses a different modification. Here, the stars are in these hierarch cal 
Jacobi coordinates, but planets are referenced purely to their primary. This permits close encounters to be implemented within each planetary 
system, but requires the small interaction Hamiltonian to be split into two separate parts. The system is evolved using leapfrog still, and the 
method is 2nd order so long as the ord er the two interaction terms is symmetric. Both symplectic integrators can be applied to the problem 
here, but that of [Chambers et alj l l2002n is favoured due to the ease of implementing close encounters. Note that, in the limit of test particles 
only or a single planet, they are identical coordinate systems. 



Table 1. Labels for the basic planetary orbital types in multiple star systems as defined in this work and compared to those o f iDvorakI Vmk . 

Orbit Triple Binary Dvorak 

description system system (1984) 
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This method requires some work to be extended to a triple system. First the coordinate system needs to be defined, as shown in Figure[2l 
S(A) planets are taken relative to star A, S(B) to star B and S(C) to star C. The barycentre of star B and its planetary system is then referred 
to the barycentre of star A and its planets (similar to lBeusg ( 120031) 's hierarchical coordinates). S(AB)-P planets are taken as relative to the 
centre of mass of the binary and all planets therein. The barycentre of the S(C) planets and their star is taken as relative to the binary and 
its circumstellar and circ umbinary planetary sys tems. Finally, P type planets are referred to the centre of mass of all the other subsystems. 
Following the notation of Chambers et alj i 2003), the transformed coordinates X are related to the inertial coordinates x by 

mtXA = niAXA + niBXB + mcxc + ''^^rnsAjXsAj + ''^^'msBjXSBj + ''^^'msCjXsCj + ''^^'msPjXsp. + mp. xp. 

j 3 J J j 
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(1) 



where subscripts A, B, C label the stars and SAi, SBi, SCi, SPi and Pi label the planets, rrit is the total mass of the system and m-b is 
the mass of the inner binary stars and their planets. To derive the full system of equ ations is somewhat inv olved, but individually (i.e. for 
one orbital type only) they are readily obtainable. For example, using the method of [Chambers et alj ( 120021) . the conjugate momenta for a 
hierarchical triple with S(AB)-P planets only are 
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(2) 



where here mb is the mass of the inner binary only. The transformed Hamiltonian is split as follows 
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where = mAruB /rrih and ~ mc{mt — mc)/mt are the reduced masses of the inner and outer binaries respectively and s — 
rrijXj/ (mt — mc) - The Hamiltonian Hkcp represents the Keplerian motion of the stellar orbits and the Keplerian motion of the planets 
about a fixed point. The recommended method by , Wisdom & Holman. (1991.) is to use the f and g functions of Danby ( 1988) to evolve the 
system under this term. 

The Hamiltonian Hint contains interaction terms dependent only on position. Hamilton's equations can be used to derive the accelera- 
tions on each object due to this term, and these can be analytically integrated to evolve the system over the interval dt. These accelerations 
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Note that in these equations V is a pseudo-velocity and equal to the transformed momenta divided by the mass of the objects they describe. 
That is, Vfl is Pfl /mb and so on. The final Hamiltonian is -ffjump and gives 



dt 



N 



j_ 

mb 



(5) 



A similar method can be used to obtain the corresponding equations for the other four orbital types, a nd these are give n in the Appendix. 
Further details of the derivation for a more gene ral system with more th an one orbital type are given in Verriei i in^prep.'). Close encounters 
can be included in the same way as described bv lChambers et al.l ( 120021) , as well as variable timesteps, but are not implemented here. 

The method outlined above was im plemented separa tely for each planetary type in a stand alone program MOIRA0, the testing of which 
is also described in the Appendix and in lVerrien din prep.!) . 



4 TRIPLE STAR SYSTEM STATISTICS 

It is helpful to determine the statistical properties o f known hierarchical triple systems, so simulations can be run that are comparable to real 
systems. The Multiple Star Catalogue jTokovinirj ri997) contains 541 reasonably complete entries for such systems, and histograms of this 
data are plotted in Figure |3] This catalogue is compiled from a variety of sources so contains a number of biases from different observing 
methods, as discussed by the author, but should be sufficient here. 

^From the histograms, it can be seen that the stellar masses tend to be near solar, with the inner binary mass ratio being generally less 
than 2 and the outer between 1 and 3, although sometimes being larger than the total mass of the binary. The semimajor axis of the inner 
binary tends to be less than 2au and the outer, although often many hundreds of astronomical units, generally concentrated around a few 
tens. There seems to be a big step in the distribution at about 15au. The ratio of the semimajor axes ranges but seems to be slightly peaked 
towards a factor of 10 or less. The inner eccentricity is peaked around zero, and a plot of the inner eccentricity as a function of semimajor 
axis shows that this does not only occur for close inner binaries. The outer eccentricity is fairly uniform from to about 0.7. A plot of the 
oute r eccentrici t y as a function of semimajor axis also shows no apparent trend. 

TokovininI (1997) considers the catalogue to be complete to a distance of lOpc. However, this includes only 14 objects with all stellar 



orbital elements known, so to determine if the data is representative it was ranked by distance and a comparison made of the first 250 entries 
to the entire sample. No major difference was apparent between the two, so the sample can be considered to give a reasonable overview of 
triple systems. 



^ The Moirai are in Greek mythology the three Fates. Given the tradition of naming planetary objects after gods and heros from classical mythology, and that 
the stars in these systems will be the largest influence on the orbital evolution of such objects, it seemed appropriate to name them (and the code) after the 
goddesses that were supposed to control the fate of men and gods alike. 
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Figure 3. The statistics of the orbits of stars in triple star systems, plotted using data from lTokovinitjfT99% . The top left panel shows a histogram of masses of 
each component (A black, B blue and C red-dashed), while the bottom left shows the mass ratios of the inner and outer binaries (mA/'mB black, mi,i„/mc 
blue). The top middle panel shows the semimajor axes of the inner (black) and outer (blue) orbits, while the bottom middle panel shows the ratio of these. The 
top left panel shows the eccentricities of the inner (black) and outer (blue) orbits, and the bottom left panel shows the ratio of the two. 




Figure 4. Details of the sets of simulations with inner binary parameters as = lau, eg = 0, = ras = lAf© and outer binary parameters ec = 0.6, 
mc = 4Mq and = 20 to 100 au. The left and centre panels show the conservation of total energy and angular momentum for these eight simulations (note 
that the change in angular momentum is a factor of 10^^ smaller than that in energy, as in symplectic integration schemes angular momentum is conserved 
to machine precision). The right hand panel shows the variation of the semimajor axis of star C throughout the simulation. The different colours indicate the 
different initial values of the semimajor axis of star C, as apparent from the right panel. 



5 NUMERICAL INTEGRATIONS: S(AB)-P ORBITS 

In this section, the numerical simulations of S(AB)-P type orbits are described and the results presented. The stability of this region between 
the inner and outer binary can be determined using grids of test particles. This is an efficient method to map out stable radii for different 
stellar orbital parameters. Since the stellar orbits and masses constitute a 7x7x7 dimensional space some simplifying assumptions are needed 
to make the investigation feasible. First, all orbits are taken as coplanar. Second, the initial orientation of the three stars is fixed, so that 
they start initially aligned all at pericentre. The effect of these assumptions will be discussed later. This leaves the masses, eccentricities 
and semimajor axes of the stars as free parameters. As the dynamics will scale with certain combinations of these, the system should be 
characterised by ratios of some of them. 
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Figure 5. The stability radii for S(AB)-P type test particles as a function of outer binary semimajor axis and eccentricity for eight different mass ratios. The 
colours indicate eccentricity of star C and the panels are labelled with the mass ratio. Solid lines show the first and last fully stable radii, and crosses show any 
unstable locations within this annulus. For comparison the location of the inner and outer critical semimajor axes predicted from lHolman & Wiegerll j 19991) 
are shown as dashed lines. 
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Figure 6. The stability radii for S(AB)-P type test pailicles as a function of outer binary semimajor axis and eccentricity, this time scaled to the semimajor 
axis of star C, for the smallest and largest mass ratios studied. The symbols are as for Fig|5] and it can be seen that the location of the outer stabihty boundary 
scales with ac ■ 




Initial semimajor axis (au) Initial semimajor oxis (au) Time (years) 

Figure 7. Some examples of test particle evolution for simulations with fic = 0.09. Different colours indicate different values of ac (see Figure|4). The 
left panels show average test particle survival time at each initial semimajor axis for = 0.0 (top) and = 0.2 (bottom), indicating that for the zero 
eccentricity case the stable region is well defined, and that that the dynamics scale with ac- The island around the 5:1MMR in the ec = 0.2 case is clearly 
visible, and the location of the resonance overlayed as a dashed line. The middle panel shows the fraction of test particles surviving at each initial semimajor 
axis and The right panel shows test particle decay rates. A fast clearing out of the unstable regions is seen. 



Given the discussion in Section |4l the stellar systems investigated is as follows. The ratio of the masses of the inner binary stars 
mA/mB is varied from 1 to 2, and the ratio of the outer binary mc/rrih from 0.1 to 2. Note that these are not the mass ratios defined in 
Holman & WiegertI Jl999 *), which are fiB ~ mB/m^j and nc ~ mc / (jnb + mc). The mass of the inner binary is fixed at 2 Mq and the 



eccentricities of both orbits range from 0.0 to 0.6. The semimajor axis of the inner binary is varied between 1 and 5 au and the semimajor 
axis of the outer binary between 20 and 100 au. The effect of changing these six parameters (^b, fJ,c, cb, ec, aB and ac) is the purpose of 
this investigation. 

For each of these stellar configurations, a grid of test particles is added. Particles in this grid are spaced evenly in radius and longitude, 
and all on initially circular orbits. A semimajor axis is defined to be stable if all particles starting there remain stable themselves for the length 
of the simulation. Individual particles are considered unstable if their orbits are not bound to the barycentre of the inner binary or if they pass 
certain radial limits, as per the definition of the orbital type in Section [2] When either of these conditions are met, the particle is removed 
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Figure 8. Test particle critical outer semimajor axes for the ac = 100 au cases: data and fits. In the left panel is the simulation data as a function of mass 
ratio, with symbols as before. Th e middle panel show s the 4 parameter fit to the data: now crosses correspond to the values to be fitted, the solid line is the fit 
and the dotted lines the results of lHolman & Wiegerj j 19991) for comparison. Note that the critical semimajor axis is taken as the innermost radii within the 
stable island not the line shown in the left panel. The right panel shows the fit as a function of eccentricity. 




Figure 9. Test particle stability as a function of initial semimajor axis, as for Figure|5] in more detail for the fic = 0.09, fic = 0.60 and fj,c = 0.67 
simulations for various eccentricities of star C. Here additional simulations have been run for values of ac between 20 and 50 au in steps of 1 au. Symbols are 
as bef ore: the solid lines are the inner and outer stable boundaries, crosses are unstable locations within these, and dotted lines the fit of lHolman & Wiegeill 
< 19991) . 



from the simulation. For S(AB)-P orbits the radial limits are the radius of the inner and outer binaries. In practice, most particles are lost as 

their orbits become unbound. 

In their study of planetary stability in binary star systems, iHolman & Wiegert ( 19991) find inner and outer critical semimajor axes 
for test particles in P and S orbits respectively. These are empirically fitted to a function of the binary's eccentricity and mass ratio (/i = 
f7t2 / {nil + 1712 ))■ For S(AB)-P orbits in a hierarchical triple system, it is reasonable to expect there to be an inner and outer critical semimajor 
axis, primarily controlled by the inner and outer binaries respectively. If the stars are separated enough to be well approximated as two 
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Figure 10. The inner critical radius as a function of the inner binaries eccentricity and mass ratio. As for Figu re[5]the solid lines show the first and last stable 
radius, crosses any unstable points between these two and the dashed lines the fit o f lHolman & Wiegerll fT999l) . 
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Figure 11. The fit to the inner critical radius as a function of t he inner binaries eccentrici ty and mass ratio. Now shown as crosses are the locations of the 
inner first stable radii, and still shown as dashed lines is the fit of lHolman & Wiegerll )l999l) . Compared to these as solid lines is the 4 parameter fit discussed 
in the text. Colours indicate mass ratio and eccentricity as for Figure [TOl 



decoupled orbits, then these critical semimajor axes should be similar in form to those found by Holman & Wiegerl 1 19991) . There are two 
obvious ways that the system differs from this simple picture. First, as the stellar orbits will vary with time, the initial orbital configuration 
may not define the maximum extent of instability regions controlled by parameters such as eccentricity. However, for the systems studied 
here it will be shown t hat this is not an impo rtant effect as the orbits do not evolve significantly. Secondly, the instability zones in binary 
systems were shown bv lMudrvk&Wijj2006h to be due to overlap between resonances. It may occur then that there are additional unstable 
regions in hierarchical triples where the resonances due to each sub-binary overlap. 

In light of this discussion, simulations can be set up to investigate the outer and inner stability boundaries separately. First, the outer 
boundary can be studied by fixing the inner binary as two I Mp, stars in a circular I au o rbit and varying star C as discussed above. The test 
particle grids in these cases are modelled on those chosen by Holman & Wiegert ( 1999h . and extend from the radius of the inner binary to 
half that of the outer binary. They are spaced evenly in steps of 0.1 au and there are eight particles at each semimajor axis. Simulation lengths 
are 1 Myr, which corresponds to at least several thousand outer binary periods. The timesteps used are of the order of a few days, giving 
relative energy conservation at the 10~^ to 10~* level. For an example of the this, see the right most panels of Figure|4l which shows the 
energy and angular momentum variation for a case with star C set to 4Mq and with an eccentricity of 0.6. Also shown in this figure is the 
variation of the semimajor axis of the outer binary, to demonstrate that even in the most extreme case the stellar orbits are not evolving. This 
is not too unexpected as the ratio between the stellar orbits is still fairly high, even in this case. (For comparison, Eggleton & Kiseleva i 1995h 
claim that the stable radius for this hierarchical triple is 16 au, and simulations show that it disintegrates at an initial separation of around 12 
au). 

The results of these simulations are given in Figure |5] which shows test particle stability for each mass ratio. It is easy to define an 
inner (and outer) first (and last) stable radius respectively as that being the first (and last) encountered where all particles are stable. This is 
plotted for each eccentricit y of star C in the figure. In addition, there are unstable radii within these bounds, and these are plotted as crosses. 
For comparison, the fits of Inolman & Wiegerll (1 1999*) are also shown as dashed lines. Note that these are the optimum values and that there 
are uncertainties given for their fit. It can be seen that the stable regions here have fairly well defined edges, and in fact match up in most 
simulations to those of lHolman & Wiegert ( I1999I) . It is expected around the inner edge to see some odd unstable radii corresponding to the 
first n : 1 mean motion resonances (MMRs) within the region. 

The inner edge of the stable region appears mostly constant as the orbit of star C is altered, as expected. The only exceptions are the 
occasional unstable location near this boundary and a general trend at small separations of the two stellar binaries for the stability zone to be 
less than that expected from the trends of the wider triple cases. It is possible to demonstrate that, apart from these exceptions, the outer edge 
scales the outer binaries semimajor axis by plotting the data scaled to this quantity, as shown in Figure|6]for two of the mass ratio cases. 

The only unclear case in these simulations is that of fic ~ 0.09 and ec ~ 0.2. Here, there is a large unstable region that occurs just 
before the last stable radius. This effect is due to a small region of stability that appears as the last stable radius. Test particles appear to be 
stable for just about the length of the simulations. Interestingly, this island appears to be centred on the 5:1 MMR with star C. Figure|7]shows 
a comparison of the evolution of test particles in this case and that with the same mass ratio but ec ~ 0.0. These show as a function of initial 
semimajor axis average test particle survival time and fraction of surviving particles. Also shown are the decay rates in each simulation. From 
these, it can be seen that the dynamics in each ac case scale with this quantity and time. It is also clear that in the ac = 20 au simulation 
the region around the 5:1 MMR is unstable with a lifetime just less than IMyr, but that as ac and the dynamical time increase the region is 
seen as stable as the test particles are not lost before the end of the simulation. In fact, the test particles within this region all have very high 
eccentricities at the end of the simulation as compared to the rest of the stable zone. In the ec ~ 0.0 graphs it can be seen that the stable 
region is very clearly defined, with very few particles surviving for any length of time outside its boundaries. 

Apart from th i s one feature for the fic = 0.09 case, the simulation lengths appear to be sufficiently long to describe stability. 
Holman & WiegertI jl999l) use an integration length of 10^ binary periods, comparable to those used here for all but the very low mass 
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Table 2. Fitted parameters for Equation|6] the outer stability edge, compared to those of iHolman & Wiegerll jl999l) . The first column shows the results of a 6 
parameter fit to the data and the middle column the results of a 4 parameter fit. 



Parameter This work This work JMtnan^^rtegert ^999) 



ai 0.477 ± 0.001 0.466 ± 0.001 0.464 ± 0.006 

a2 -0.412 ± 0.002 -0.392 ± 0.001 -0.380 ± 0.010 

as -0.708 ± 0.006 -0.542 ± 0.002 -0.631 ± 0.034 

04 0.794 ± 0.012 0.494 ± 0.004 0.586 ± 0.061 

as 0.276 ± 0.009 - 0.150 ± 0.041 

ae -0.500 ± 0.020 - -0.198 ± 0.074 



Table 3. Fitted parameters for Equation[8] the inner stability edge, compared to those o f lHohnan & Wiegerll jl999h . The first column shows the results of a 7 
parameter fit,and the second column the results of a 4 parameter fit. 



Parameter This work This work liolman^Wiegert ^999) 



ai 


3.45 


± 


1.10 


2.92 


± 


0.12 


1.60 


± 


0.04 


02 


9.94 


± 


1.81 


4.21 




0.24 


5.10 




0.05 


as 


-6.95 


± 


1.47 


-2.67 




0.38 


-2.22 


± 


0.11 


a4 


-5.43 


± 


5.27 


-1.55 




0.29 


4.12 


± 


0.09 


05 


-14.09 


± 


4.42 








-4.27 




0.17 


ae 


6.22 


± 


6.26 








-5.09 


± 


0.11 


ar 


25.46 


± 


8.51 








4.61 


± 


0.36 



ratio and large ac cases. They also mention that the stability boundaries do not change much after a few hundred binary periods. This would 
indicate that even for the long binary period cases, the results are unaffected by this choice, and this is supported by the scaling seen with 
semimajor axis in the location of the outer boundary. As shown in Figure|7]and discussed above, the edges to the stable regions are distinct 
in most cases, and the test particle decay rates indicate that a stable situation has been reached long before the end of the simulations. Further 
support that the simulations times are sufficient is given by running the ~ 0.09 and ec ~ 0.6 cases for 2 Myr. The results from these are 
almost identical to the original 1 Myr simulations, with only a few additional unstable points appearing. 

Since the outer boundary is well modelled as a function of mass ratio and eccentricity and scales with ac for all but the smallest 
separations the ac ~ 100 au data can be plotted and fitted. This semimajor axis is chosen as it has the most detailed determination of the 
critical radii due to the larger number of test particles used. To provide a better fit additional simulations were run for six more mass ratios, 
and the data to be fitted is shown in the left hand panel of Figure[8l I nstead of using the last stab le radius, the stable locations just within any 
odd unstable radii around the boundary are taken instead. Following iHolman & WiegertI il999h , we fit this boundary with a function of the 
form 

= tti + a2Hc + a^ec + a^^c^c + "see + aefJ-c^c (6) 

ac 

Performing this fit gives th e parameters as shown in the left hand column of Table |2] Also given for comparison are the values suggested by 



Holman & Wiegert 1 1999h . This fit has a of about 2600 and is in reasonable agreement with their values. Note however that here a smaller 



range of mass ratios has been fitted, as it is not realistic to increase the mass of star C to the point where fic = 0.9. A four parameter fit to the 
first four terms in Equation|6]gives a simpler fit, also shown in the table. The value of this fit is higher at 3700 but still reasonably models 
the data, and this fit is shown in the middle and right panels of Figure[8] Note that the value is calculated using the standard formula 

(7) 

where yi is the fitted value of the function, yotsi the observed value and (Ji the error this value. Since the grid size here is 0. 1 au the stability 
boundary cannot be located to any greater accuracy and hence this was assigned as the uncertainty in each measurement. The large x^ values 
obtained for the fits reflect the complex nature of a boundary that is not easily fitted with such a simple function. However, as can be seen 
from the figures, it is still a useful approximation for quickly determining the stability properties of a system. 

As mentioned above, the inner stability boundary seems to be unaffected by changes in the orbit of the outer star except for at small 
separations of the two binary orbits. In this case, the stable region is smaller than expected and there are also far more isolated unstable radii, 
especially as the mass of star C increases. Figure |9] shows effect in more detail for the cases of /ic = 0.09, 0.60 and 0.67. Here several 
more semimajor axes for the outer star have been studied for various eccentricities. Although there seems to be a linear overlay of the two 
different boundaries there does seem to be a slight effect around the point where they overlap, becoming more pronounced as the mass ratio 
and eccentricity increases. There also seems to be a definite unstable island in the centre of the region, most noticeable in the ec ~ 0.4 cases. 
This moves outwards as the position of the outer star is increased, and must be a resonant feature. Note that if the mass of star C is set to zero 
for the highest mass ratio and eccentricity case then the stable region has a clearly defined inner radius at 2.3 au and there are no unstable 
particles beyond this boundary, indicating that the structure seen in these graphs is a direct result of the combined effects of all three stars. 
Whether this is due to resonance overlap is unknown. 
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Figure 12. The stability radii for S(AB)-P type test particles as a function of outer binary semimajor axis and mass ratio for different eccentricities of both 
stellar orbits. The colours indicate eccentricity of star B and the panels are labelled with the eccentricity of star C. Solid lines show the first and last fully stable 
radii, and crosses show any unstable locations within this annulus. 



A more detailed investigation of the inner boundary can be carried out in a similar manner to the outer edge, by fixing star C at 1 Mq and 
in a circular 50 au orbit and varying the inner binaries mass ratio and eccentricity, as discussed earlier. The radius of the inner pair was kept at 
1 au and their total mass as 2 Mq . The eccentricity of this binary was varied from 0.0 to 0.6 in steps of 0.2 again, and the mass ratio mA /ruB 
varied from 1.0 to 2.0 in steps of 0.1. Figure [To] shows the results from these simulations, with the locations of t he first and last stab l e radi i 
plotted as functions of the inner binaries eccentricity es and mass ratio fiB- Also plotted again is the fit given bv lHolman & Wiegerll Il999l) 
to the critical radius for P type planets in binary systems. As expected, the outer stability boundary seems unaffected by changes to the inner 
binary pair. 

There are some unsta ble points between the two stability boundaries, most notably around 3.2 au for the simulations with Cb ~ 0.0. 
Holman & Wiegertlll999l) find unstable islands appearing at the first n : 1 beyond the critical semimajor axis in this configuration. However, 



the location of these unstable radii here is well beyond the first n : 1 MMR in the stable region. In addition, running the hb = 0.34 and 
es ~ 0.0 case without star C reveals that these locations are now stable. This would indicate that this is again an effect due to the combination 
of all three stars. Because of this the inner edges location is simply taken as the same as the first stable radius. 

The location of this boundary app ears to be a very weak fimc tion of mass ratio and an approximately linear function of eccentricity, and 
also agrees well with the predictions of Holman & Wiegerl ( 19991) . For this boundary, they fit a function of the form 

'^in. , , 2 , , .2,22 /o\ 

= ai + a2es + aae^ + 04/1^ + a^eB^JiB + asMs + aresfJ-B (o) 

as 

where the constants are given in the third column of Table [3] The terms up to fourth order are included to fit a variation in the position of 
the critical radius at smaller mass ratios than those considered here. Fitting this function to the results here does not produce well determined 
coefficients despite a low value of about 29, as shown in the first column of the table. In fact, a better solution is obtained by only including 
the first four terms, as shown in the middle column of the table. The value for this fit is slightly higher at about 42. This second fit is 
plotted in Figure fTTl and seems to describe the data well, although it seems to slightly overestimate the size of the stable region. Despite the 
smaller values here the fitted parameters are less well determined than those for the outer edge. This is a reflection on the relative size of 
the test particle grid spacing compared to the size of the inner binaries orbit. Here the boundary is determined to within a tenth of the relative 
separation of the stars, while for the outer edge it was determined to within a thousandth of the size of the binary's orbit. 

The parameter space investigated so far is somewhat limited. The stars orbital longitudes have been ignored, assumed to be a minor 
influence on stability, the sets of simulations have always kept one star on a circular orbit, and all objects have been taken as coplanar. Brief 
investigations of these three extensions to the parameter space can be made. 

Firstly, the assumption that the initial longitudes of the stellar orbits has little effect on the stability boundaries was tested by running 
the inner edge simulations with mass ratios fiB =0.33 to 0.40 with a different initial longitude of the inner binary pair. The results from these 
simulations were almost identical to those of the initial set, providing some evidence that this assumption is valid. 

Next, the effect of both stars having non-circular orbits was investigated. The semimajor axes of the inner and outer binary were kept 
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Figure 13. The stability radii for S(AB)-P type test particles as a function of outer binary semimajor axis and mass ratio for different eccentricities and 
inclinations of star C. The colours indicate inclination and the panels are labelled with the eccentricity. Solid lines show the first and last fully stable radii, and 
crosses show any unstable locations within this annulus. 



at 1 and 50 au, and the eccentricities of both varied in steps of 0.2 from 0.0 to 0.6. The mass ratio of the outer binary only was varied as 
before. The results of these simulations are shown in Figure [12] Each panel shows the location of the stability radii and unstable points as 
for all values of es and one given value of ec- There is almost no change in the position of the outer stability boundary, but as both stellar 
eccentricities and the mass ratio increase, the inner boundary starts to move outwards. The stellar eccentricities are still not varying to any 
significant extent and the additional instability must be due to the combined perturbations of all three stars. 

Lastly, the effect of inclination was studied. The semimajor axes were kept as 1 and 50 au, bb was set to 0, and the outer mass ratio 
and eccentricity varied as before. Sets of simulations where then run for inclinations of the outer star of 10°, 20° and 30°. The test particles 
were kept coplanar with the inner binary. This is a rather limited investigation, as any dependence on the longitude of ascending node has 
been ignored, and only a small range of inclinations included. However, higher inclinations will be subject to the Kozai instability, causing 
large variations in the stellar orbits, which is expected to rapidly destabilise test particles. In these simulations the stellar mutual inclination 
remains fairly constant and their orbits are stable. Figure [13] shows the stability boundaries for these simulations, each panel comparing the 
different inclination results for a different value of ec- There is little change in the inner boundary but the outer edge moves somewhat. 
Interestingly, for the ec = case higher inclinations are more stable. If the test particles are st arted instead coplana r with the outer binary 
the stability is similar, although not identical. These results are consistent with the conclusions of lPilat-Lohinger et al.l ( 120031) , who show that 
inclination is not a significant effect on the stability of P type planets in binary system. 



6 CONCLUSION 

The main achievement of this paper is the formulati on of a symplect i c inte grator algorithm suitable for hierarchical triple systems. This 



extends the algorithm for binary systems presented by [Chambers et alj ( 120020 . The positions of the stars are followed in hierarchical Jacobi 



coordinates, whilst the planets are referenced purely to their primary. Each of the five distinct cases, namely circumtriple orbits, circumbinary 
orbits and circumstellar orbits around each of the stars in the hierarchical triple, requires a different splitting of the Hamiltonian and hence 
a different formulation of the symplectic integration algorithm. Here, we have given the mathematical details for each of the five cases, and 
presented a working code that implements the algorithm. 

As an application, a survey of the stability zones for circumbinary planets in hierarchical triples is presented. Here, the planet orbits an 
inner binary, with a more distant companion star completing the stellar triple. Using a set of numerical simulations, we found the extent of 
the stable zone which can support long-lived planetary orbits and provided fits to the inner and outer edges. The effect of low inclination 
on this boundary is minimal. A reasonable first approximation to a behaviour of a hierarchical triple is to regard it as a superposition of the 
dynamics of the inner binary and a pseudo-binary consisting of the outer star and a point mass approximation to the inner binary. If it is 
considered as two decoupled binary systems, then the earlier work of Holman & Wiegert (1999) on binaries is applicable to triples, except in 
the cases of high eccentricities and close or massive stars. 
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Figure 14. Widths of the circumbinary stabihty zones for known triple systems in au and as a percentage of the area between the two sub-binaries, calculated 
using the four parameter fits derived in Section|5] Note that using the criteria of lHolman & Wiegerll (1993) gives almost identical results. 



The impli c ation here is that the addition of a stable third star does not distort the original binary stability boundaries. As mentioned, 
Mudrvk & j2006h have shown that overlapping sub-resonances are the cause of the boundary in the binary case. It is reasonable to 



expect that in triples the same process is responsible, and the similarities between the binary and triple results support this theory. It is also 
expected that there is a regime in which the resonances from each sub-binary start to overlap as well, further destabilising the test particles. 
Evidence of this is the deviation from the binary results when the stars are close, massive and very eccentric, when resonances would be both 
stronger and wider. Since the parameter space investigated was chosen to reflect the observed systems it would seem to be a reflection on the 
characteristics of known triple stars that most lie in the decoupled regime. The relatively constant nature of the stellar orbits in the simulations 
is however a consequence of the test particle orbits being destabilised long before the stars are close enough to interact. By extension of all 
these arguments, it is expected that the binary criteria can be used to fairly accurately predict the stability zones in any hierarchical stellar 
system, no matter the number of stars. 

The r e sults p resented here can be used to estimate the number of known hierarchical triple systems that could harbour S(AB)-P planets. 



Ine r e sults p resented tiere can be used to estimate trie number ot known hierarchical triple systems that could harbour i(Alj)-F planets. 
Tokovinii] 1997h lists 54 systems with semimajor axis, eccentricity and masses for both the inner and outer components. The mutual 



inclinations of most are not well known, but there are nine systems listed for which this angle can be determined. For five of these it is less 
than 15°, two are around 40° and two are retrograde. A lthough a small sample it suggests that there are systems that fall within the low 
inclination regime investigated here. Using the criteria of Holman & Wiegert ( 19991 ) and those found here for the position of the inner and 
outer critical semimajor axis the size of the coplanar stable region for each of these triples can be calculated. This can be considered an 
upper limit, since it is likely that very non-coplanar systems and those with significant eccentricities for both binary components will further 
desta bilise planets. Of the 54 sys tems 13 are completely unstable to circumbinary planets according to the four parameter fits (compared to 1 1 
using^Holma n & Wiegert 1 1999V s criteria). Figure [14] shows a plot of the width of the stable region for the remaining systems. Interestingly, 
the majority seem to have very small stable zones, with 16 smaller than an au. Whether this is a feature of triple systems, or an observation 
bias is not apparent. It does indicate though that circumbinary planets are unlikely to exist in at least 50 % of observable systems. 
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APPENDIX A: APPENDIX 

In this appendix are details of the symplectic integration method derived in Section[3]for the orbital types P, S(A), S(B) and S(C). The split 
Hamiltonian in each case is presented, and the evolution under these is readily obtainable using the method given in Section[3] Also presented 
are details of the testing of the implementation of this symplectic integrator. Note that throughout this appendix mt is the total system mass 
and mb the total inner binary mass, including planets where relevant. 
For S(A) orbits, the Hamiltonian is 

H — Hkcp + ^^Int + J^Jump 

JV 



-HKcp ~ r 5 5 h 




Rb Rc ^ \ Irni Ri 



Rb \Xb+s\ 



niB mcm-b 



\Xb — Xi + s\ \mbXi — jubXb — rribXc — rnbs\ I 

1 niA ruB \ 

Rc \mBX B + rrihXc + mhs\ \mhX c — rn(^A+p)X b\ j 

H]ump ~ — — (Al) 
2mA 

where m^A+p) is the mass of star A and its planets, fitin is the reduced mass of the inner binary, including the planets, jj,tot is the reduced 
mass of the outer binary and s — rrijXj /m(A+p). 
For S(B) type planets only, the equations are 



H = Hkcp + i/int + i/,7u 



mp 



^ _ Pb _^ Pc _ Gfj.bi„{mt - mg) ^ G^totmx ^ / P^ _ GmBmj 
2^bi„ 2fitot Rb Rc ^ \ 2mi Ri 



Hi, 



N 

Gmiirij I m^B+p) 




mc 



Xi + qXb — Xc — s\ 



Hjump ~ — (A2) 

2ms 

where m(^B+p) is the mass of star B and its planets, is the reduced mass of the inner binary, including the planets and fitot is the reduced 
mass of the outer binary. Here s = X/ ^i-^ i /''^(B+p) now, and q = mAjijnx ~ mc) and r — m^B+p)/ {mt — mc)- 
The corresponding Hamiltonian for S(C) type planets is 

H — Hkcp + Hint + Hjump 

P% Pa GmAmB Gm^c+p)mb / Gmcmi 

ijKcp = h h 



2^bin 2iitri Rb Rc ^-^ \ 2mi Ri 
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j = l j>i 
N 



Ri 



Rc 



- Xc + rXi 



Gmc 



+ 



\Xc + rXB-s\ 



TUB 



+ 



Xc 



■ qXs — s\ 



\Xi + Xc — qXb — s\ 



H 



Jump 



2mc 



(A3) 



where Ubin is the reduced mass of the inner binary stars and fitri is the reduced mass of the outer binary, including the planets. Here now 
s = ^ nijX j /m^c+p)^ Q = rnA/rrih, r — niB /nih and mf^c+p) is the mass of star C and its planets. 
Finally, for P orbits, the equations are 



H 



Hkcp + Hint + Hjun 



Kcp 



p2 



2^6 



+ 



Gflbinmh Gj-Ltrimtri 



2flt 



Gmirrii 



^ — ' ^ — ' Rij 

j = l j>i 

N 

— Grrii, ^ ^ nii 

N 

+ Gmtri ^ "1 

\EP. 



Rb 



+ Gmcnib 



Rc 



+ 



N 

E 



p2 

2mi 



Gmtrinii 
R. 



TUA 



rriB 



IniBXB + mhXc 



IruAXB - mhXc 



mA 



niB 



|mbXi + mflXs + fitriXcl \mi,Xi — ttiaXb + /itri^cj 



mc 



mtr 



m^Xc 



(A4) 



where Ubin is as before, utri is the reduced mass of the outer binary's orbit and mtri is the mass of all three stars. 

These differing orbital cases are all dealt with separately at present. To test the method and its implementation comparisons of the 
simulations of systems of various levels of complexity were made with the literature. As the equations given above allow mc to be set to 
zero, the conser vation of the Jacobi con stant of test particles in the Circular Restricted Three Body Problem was calculated and compared 
to that given by [Chambers et alj ( 12003) for their binary syste ms. The evolu tion of various triple star system was compared to studies of 
hierarchical tr iples in the li terature, for example those given in lBeuslI ( 120030 . Simulations of discs of test particles were also compared to 
those given bv lBeuslllioOsI) for circ umbinary discs. M ore complex systems of planets and test particles were also compared to the results of 
a standard Bulirsch-Stoer integrator dPress et al ]|2002h . In all cases excellent agreement was found. 
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